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ABSTRACT 


This report covers work on two unrelated tasks. The 
first was an investigation of the formulation of the equations 
for non - uniform unsteady flows, by perturbation of an 
irrotational flow to obtain the linear Green's equation. 

The resulting integral equation was found to contain a kernel 
which could be expressed as the solution of the adjoint flow 
equation, a linear equation for small perturbations, but with 
non-constant coefficients determined by the steady flow 
conditions. For the uniform flow case, this kernel was 
found to limit to the doublet form commonly used in formulating 
the flutter problem. It is believed that the non-uniform 
flow effects may prove important in transonic flutter, and 
that in such cases, the use of doublet type solutions of the 
wave equation would then prove to be erroneous. The second 
task covered an initial investigation into the use of the 
Monte Carlo method for solution of acoustical field problems. 
Computed results are given for a rectangular room problem, 
and for a problem involving a circular duct with a source 
located at the closed end. In both cases, results appear 
to be within statistical expectations, although statistical 
deviations were large because computer time was limited. The 
most severe limitation to further applications of the method 
is the need for a suitable method to handle acoustical 
diffraction, as in a shadow zone, or near a window in a room. 
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INTRODUCTION 


Work under this grant started on October 1, 1969 and 
was terminated on September 30 , 1973. The main effort was 
concentrated on two tasks: (A) Numerical Methods of 

Solution of Non-Linear Unsteady Transonic Flows, undertaken 
entirely by the author, and (D) Application of Monte Carlo 
Methods to Acoustical Analysis, undertaken as a doctoral 
dissertation by Balakrishna Thanedar, with the author acting 
as adviser. Two other tasks were briefly worked on, but were 
abandoned before any appreciable effort had been expended. 

These were: (B) Dynamic Interaction of Tracked Air Cushion 

Vehicle and Guideway; and (C) Computer Model for Aircraft 
Ride Environment Analysis. 

Work under task (A) became divided into two major subtasks. 
The first was continued up to the termination of the project, 
and consisted of investigations into methods of formulating 
the nonlinear transonic flow problem. It has resulted in 
a Note 1 to the AIAA Journal, "The Integral Equation for Small 
Perturbations of Irrotational Flows", in which it is shown 
that the Green's Function required for the formulation of the 
integral equation relating normal flow velocity to local 
velocity potential is the unit solution of the adjoint wave 
equation. Hitherto, it has generally been assumed that it is 
the unit solution of the wave equation. 
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The second subtask concerned the formulation of the 

flutter equation in terms of velocity potentials. It was 

argued that the nonlinear transonic problem might be approached 

more readily in terms of the velocity potential. However, 

before this formulation could be used with confidence in the 

nonlinear solution, it would be necessary to check it against 

known results for the linear case. A computer program was 

prepared, using not only the "downwash-velocity potential 

method" (often called the "integrated potential method") 

but also introducing the concept of "aerodynamic elements". 

This effort was transferred to NASA Grant No. NGL-47-005-098 

"Numerical Methodology for Flutter Analysis and Optimization 

of Aircraft Structures" during 1971. It has since resulted 
2 

in a paper published in the AIAA Journal; "Downwash-Velocity 
Potential Method for Oscillating Surfaces". Two more 
papers are planned, and a final report on the work under this 
grant is in review. The present status of the method is 
that it has been developed for out-of-plane polygonal 
elements in subsonic and supersonic flow. 

Work under task (D) had its first results in the 
demonstration of a Monte Carlo solution for a rectangular 
room. This was presented to the April 1971 meeting of the 
Acoustical Society, and will be published 3 shortly in their 
Journal as "Monte Carlo Applications to Acoustical Field 
Solutions". Since then, a solution for a circular duct was 
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obtained. Both of these solutions were described in a 

. . 4 

doctoral dissertation by Thanedar , "Monte Carlo Investi- 
gation of Transient Acoustic Fields in Partially or Completely 
Bounded Media". 

Discussions of the two tasks are given separately in 
Parts I and II of this report. 



PART I 

METH0DS 0F SOLUTION OF 
NONLINEAR UNSTEADY TRANSONIC FLOWS 

DISCUSSION 

The aerodynamic disturbances which are transmitted from 
pornt to point on an oscillating wing depend, for their 
velocity, on the local velocity of sound and on the local 
airflow velocity, it might therefore appear obvious that 
local variations of the speed of sound would result in 
appreciable modification of the aerodynamic behavior of the 
wing when the airflow is locally sonic, interest in this 
problem has been spurred by poor agreement between experiment 
and theory in transonic flows as compared to moderate to good 
agreement in subsonic and supersonic flows. The nonlinear 
transonic problem has been studied by a number of authors, 
generally as a linear perturbation. 

Two major considerations spurred the work performed 
under this grant, , 1 , that the generally accepted linearised 
form of the integral equation was open to suspicion, as it 
seemed to imply that a pressure disturbance on a wing anticipates 
the local normal flow resulting from a deflection of the 

surface; and (2) that the downwash velocity potential formulation 
IS a potential building block of any numerical formulation, 
because the kernel of the integral equation directly relates 
velocity potential at one point to normal flow velocity 


4 
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(downwash) at another. In the downwash-pressure formulation, 
on the other hand, the kernel relates an integral of the wake 
to infinity behind one point to the normal flow velocity at 
another. 

A tentative approach to the development of a transonic 
method was formulated as follows: 

(a) Develop a subsonic uniform flow method based on the 
downwash- velocity potential formulation. 

(b) Develop the same for the supersonic case. 

(c) Develop a ray tracing method to permit calculation 
of the time lag and attenuation between pairs of points. 

(d) Develop a steady flow solution for finite thickness 
wings. This is needed to provide a medium through which 
rays can be traced. 

(e) Combine into one program. This would perform the 
following: 

(i) Compute local steady state flow velocities and 
Mach numbers 

(ii) Perform ray tracing calculations connecting pairs 
of points 

(iii) Form the aerodynamic matrix by substituting 
actual attenuation and time lag values, calculated by ray 
tracing, into kernels of expressions developed from linear 
equations . 
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(iv) Use the aerodynamic matrix as appropriate to 
complete the solution of the problem. 

Steps (a) and (b) were initiated under this grant and 
have since been achieved in separate efforts, as noted in 
the Introduction. Step (c) has been taken to the point that 
the correct integral equation has been formulated. The 
next move requires ray tracing with the adjoint equation, 
which has not yet been attempted. Steps (d) and (e) have not 
been initiated. 

Different methods of analysis which have been proposed 
are discussed in the following. section, and work on the formu 
lation of the integral equation, as required for step (c) , 
is described in Section 1.3. 
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1.2. REVIEW OF METHODS OF ANALYSIS 

Ray Tracing 

As will be shown in the following section, the correct 
formulation for the flow velocity in the normal direction 
is given by Eq. (1-53) . This expression includes the Green's 
function Gr , which must satisfy the adjoint wave equation, 
as expressed by Eq. (1-30) . However, it has often been 
assumed that the Green's Function and the source solution 
defined by Eq. (1-43) are identical. The precise effect of 
this assumption is not understood at this time, but, 
as shown in the next section, they are identical in the 
uniform flow case, and have the following values 


& - tKT*)<3R. + u(t a ;g a 

where the subscripts R, A stand for retarded and 
advanced, respectively, and is a cutoff function which 

equals unity when X is real and positive, zero otherwise. 


Also 


Gr, = - e LtOT V 4- FT ft 
Ga = - e." ^ Ta / ZfrrR 


where R - J CX p 

The points x, y, z and ^ are the 'receiving' 

point iff" and the 'disturbing point’ jp , respectively. 

The terms , Tft are the times taken by the retarded 

and advanced waves, respectively, to go from the disturbing 
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point to the receiving point, and are given by 

Tr. ' M Cn-J) } /pc. 

T A = { n('t-'f) } />\ 

where ^ - I — / v '| 
and c = speed of sound. 

It is apparent from the above expressions that only 
the retarded wave arrives at the receiving point in a subsonic 
flow, and that both waves arrive in supersonic flow only if 
the receiving point is located in the aft Mach cone from the 
disturbing point. 

An obvious approach to a solution of the nonlinear 
problem would be to determine values for ^ t R , h*) 

and c appropriate to the actual non-uniform solution. 
Evaluation of these quantities would require some form of 
solution of the non-uniform wave equation, possibly by 
ray tracing. 

Several possible ray tracing approaches have been 
considered, but since it has been determined that rays should 
be traced with the adjoint equation, as demonstrated in the 
next section, no firm understanding has been developed of 
how this ray tracing is to be done. 

The approach described here was originally suggested 
by Andrew and Stenton , the idea of determining and 
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and substituting new values for *7^ and '7^ having been 
suggested by Landahl 6 . 

Small Perturbation Methods 

A summary by Bland 7 shows that five different approaches 
to the transonic aerodynamics problem are based on the small 
disturbance potential equation given by Landahl 9 . Assuming 
that the velocity potential can be expressed as the sum of 
a steady and an oscillating part, as in Eq. (1-14) , a linear 
differential equation is derived with variable coefficients. 

This is also true of Eq. (1-20) , but the two equations are 
different, presumably due to approximations introduced into 
Landahl' s, which is valid only for transonic flows. 

Finite-Difference Method 

This method has been investigated by Ehlers 9 for the 
two-dimensional case, with promising results. Within the 
limitations of computer storage, it should be possible to 
obtain results by this method, regardless of how the perturbation 
equation is formulated. 

Local Linearization Procedure 

This analytical method has been investigated by Stahara 
and Spreiter 10 , and has been developed from methods used for 
the solution of steady problems. 

Layered Medium Analysis 

This has been investigated by Revell 11 . The region 
around the wing is divided into subsonic and supersonic areas, 
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which are further subdivided into horizontal layers in which 
local speeds of sound and Mach numbers are assumed to be 
constant. Linear solutions in each layer are matched at the 
boundaries to form complete solutions. Numerical results 
are not yet available. 

Lifting Surface Element Method 

This has been investigated by Cunningham 12 . The 

surfaces are divided into smaller elements, each assumed to 

be in uniform flow, but with a flow condition different to 

tha adjacent elements. The normal velocity components at 

collocation points are then calculated much as in many of 

the variants of Watkins ' 13 kernel function method, more 

especially as in the collocation method developed by 
14 

Cunningham , except that the local flow conditions are used 
in each case. 

Modified Sonic Box Method 

This method, investigated by Ruo, Yates, and Theisen 15 

is a modification of the sonic box method of Rodemich and 
1 6 

Andrew , in which a transformation is used to replace a 
nonconstant Mach number by a constant one. A flutter calculation 
on a delta wing using this method showed the correct trend 
as thickness was varied. 
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1-3 THE INTEGRAL EQUATION FOR TRANSONIC FLOWS 

Introduction 

The correct form of solution of the problem of oscil- 
lating airfoils in non-uniform irrotational flow is examined 
by linear perturbation. 

It is shown that the kernel of the integral equation 
giving the unknown velocity potential perturbation on the 
surface, as a function of the known normal flow or "downwash" 
perturbation, is the second derivative of a Green's function. 
It is also shown that an approximate solution of this function 
might be sought by reverse ray tracing with the adjoint wave 
equation, i.e., by starting at the "downwash" point, or 
"receiving" point, and proceeding towards the "potential 
perturbation" point or "disturbing" point. 

It might appear, intuitively, that the kernel could be 
obtained from a source solution, in which a ray is traced 
from the "disturbing point" to the "receiving" point. This 

is all the more plausible as it corresponds to the method 

. . 13 

used m the literature to obtain the kernel in the uniform 
flow case. However, in the latter case, either method would 
lead to the same result because the wave operator is 
Hermetian (i.e. its transpose is its adjoint so that its 
adjoint is also the operator for a reverse flow) . 

Unfortunately, the non-uniform flow wave operator is not 
Hermetian, so that the direct ray tracing technique leads to 
incorrect results for non-uniform flows. 
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Coordinates 

kf - Coordinate Vector of "receiving” or "downwash" point 

Jp = Coordinate Vector of "sending", "potential perturbation" 

or "dummy" point 

Y\. = Normal coordinate to bounding surface at point /f * 1 

La* = Normal coordinate to bounding surface at point 

V = Normal coordinate to one-sided surface 3 at y 

- Coordinate parallel to the flow 
^ = Dummy for 

Symbols 

C = local speed of sound 
Cy = Green's function 
P = Pressure 

p = amplitude of pressure perturbation 
= Flow velocity vector 
S = One-sided surface 
■*'£) — Bounding surface 
t - time 

"V - Steady flow velocity vector 

■O' = Volume bounded by S 

= Velocity potential and steady component 

= amplitude of velocity potential perturbation 

A 

jjjt, = Unit normal to surface 
$> — Density and steady component 

Cl) - Radial frequency 

p - Amplitude of pressure potential perturbation 
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Operators 


= Wave operator and adjoint 
TP = Bilinear concommitant 
V = Del operator 

A = Difference across a discontinuity 
Subscripts 

= Indicate variable of differentiation 
A = Airfoil surface, as opposed to wake 
L = Lower 
IF = Upper 

Potential Equation for Irrotational Adiabatic Flow 

For an irrotational flow, the velocity vector Q can be 


derived from a potential by 

« in ; 

then V ^ - O (1^2! 

Because the flow is adiabatic, the pressure is a function 
of density Jp only, so that 

P = P<Cy3 an: 

and dF/dp= cftFJ) (I _ 4; 

Writing the continuity equation as 

3^/dt 4 - o (i-5] 

and rearranging after substituting from Eq. (4) 

\ ^P 




(1-6) 
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Euler 1 s equation can be written as 


= o 


(1-7) 


From a theorem on vectors 


= ?V<f- 


( 1 - 8 ) 


Thus, substituting in Euler's equation, and noting that the last 
term disappears according to Eq. (2) , the following form of 
Euler's equation results 

tf + 


But, substituting from Eq. (1) , this can be written as 

ni - -° 


= O 


(i-io) 


which can be readily integrated to give a form of Bernoulli's 


equation 


+ a 


$ - f(H 

3 


(1-11) 


Now, preraultiply Eq. (9) by OH* , add d/dtr of Eq. (11) , and 


subtract Eq. (6) to obtain 

*§ + *£ 4 

at 2 - 3tr 






(1^12) 


With a slight rearrangement, and a substitution from 
Eq. (1) , the form given by Garrick^ results 


is + ti . 

at* at 1 


z ^ 


- c V § 


(1-13) 



15 


Perturbation Equation 

The flow is assumed to consist of a steady component 
with a small harmonic perturbation, so that 

$=§o+<t>e lwt ( i-i4, 

and 

= V + yc|) d-15) 

where both ^ and refer to the Steady flow component, and 

V " V (1-16) 

so that ^ % V “ O (1-17) 

The following expressions result from perturbing individual 
terms in Eq.(13), and dropping second order small quantities. 

&4>/dfc s - 

= V-VV+ 2iy-v)$> e Wt + [v/j-V^e^ 

V 4 |> = W'V + V Z e uufc 

Then these are substituted into Eq. (13) , two equations 
result; a steady state equation 

4 V- yv = c v- v 


(1-18) 



16 


and a perturbation equation 


/£<£>} 


(1-19) 


where 




'■* (I - 20 > 


The square bracket is intended to signify that the 
influence of the operator contained in it is restricted. 

Thus [Wjis a non-operating vector, while (V*V) i s a 
scalar operator. 

Scheme for Obtaining Green's Equation 

Equation (19) can be expressed in the more precise form 


Jl*- £<RrJ)} =o 


( 1 - 21 ) 


where IP is the vector coordinate of the "receiving" point. 

The subscript on the operator**^ indicates that its differential 
operator components, for example Vyp » operate with respect 
to fP • Next, a function (jr is assumed, satisfying the 


equation 


Xw~ r)] = S(ir-y) (i- 

The operator is the adjoint of £ , and jj£) is the 


( 1 - 22 ) 


vector coordinate of the "sending" or dummy point. 
The following expression is now derived 

Gfib} = 4>/£&l +V-P|G;4)] 


(1-23) 
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where |r is the "bilinear concommitant" . This is next 
integrated over the complete space >0“ with dummy coordinate 
vector so that 

On substitution from Eqs. (21) and (22), and on 
divergence theorem 

4 > 0r)=-J 6 jcU 

-A 

where is the unit outward normal on the surface 

bounding the space • 

The function (jr is variously referred to in the literature 
as the free-field Green's function, the unit solution, or 
when the operator is hyperbolic, the Reimann function. It 
will be referred to here simply as the Green's function, 
although this should more properly be reserved for the function 
which satisfies not only Eq. (22) , but also the boundary 
conditions of the problem. 


(1-24) 


using the 


(1-25) 
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Derivation of Green's Equation 

In order to express Eq. (23) in terms of the expression 
for given in Eq. (20) , the laws of product differentiation 
must be applied, term by term. First, these terms are con- 
sidered separately, as follows 

2 . Z 

GV(fc - V* GV(j} - CVGj-¥(i)-\V*GVct)-V'<i)VG+^VG 

V-^.G4>-cbV- %. & 

§. Lv-vf 4) - g^-' v(v-v<t>h' V-^g(v-v<|>) - b/*v4Qv~ &■ 

z 

S-Cv = v- (W] ^ - 4> V. (W] 

Then using the above results, Eq, (23) becomes 

SiY l + »-&>VV. 

= §[-+ 2uoV- ^ - ( VAY* V-2^ [Wv3]g 

+v-|g v4>-4>vg - ^ G(v-v4>) 

^ V(KV^)^-2^f>vJ 2 G4>j 

C— ' 


(1-26) 
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If the surface bounds the flow, and is nonporous, then. 


on this surface 
A 


IP-- V = o 


and 


A 

l\x- 


= 


(1-27) 


(1-28) 


where is the local normal coordinate. Thus, identifying 
the bilinear concommitant in Eq.{26), and substituting into 


Eg. (25), using Eqs.(27) and (28) 




(1-29) 


The Green's function [fep) satisfies Eq. (22) in the 


Jp space , i . e . 


fpi&CfcjO } = 


(1-30) 


where from Eq. (26) 

•v z -X 


£ - V + f 2'u»r- J-(VAV)^+V'^[VV J 


(1-31) 


Normal Flow at a Boundary 

Equation (29) gives the perturbation velocity potential 
at any point in the volume Aj~ . Generally, the problem is 
to solve for the velocity potential at any point on the 
boundary given the normal flow velocities ^ ^ At • 
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where is the local normal coordinate to the boundary 

at point IP • Thus, from Eq. (29), an integral equation is 



Assuming no radiation from infinity, the surface >£> 
consists of immersed bodies and any wake which sustains 
potential jumps. 

Uniform Flow Over Thin Airfoil 

A thin airfoil and its wake are shown iiniFigure 1. It is 
assumed to be a small angle of attack so that C and V are 
uniform everywhere, with V parallel to the %- axis, and of 
magnitude \f . 


V 



It is clear that values of on opposite sides of the 

airfoil are equal, while the local normals are equal and 
opposite. Thus the integrand containing S(|)/3jLi.cancels out, 
while, if the area of integration is transformed to a surface 
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on one side of the airfoil and wake, as shown in 
Figure 2, in the second integrand is replaced by^i(j^ , 

the velocity potential difference. If V is the normal 
coordinate, as shown, then is positive when <|> 

decreases in the positive V direction. 

V 


Suu-^t>ce -S 


r 




+ 

v> 


- c|)(_^ a 4^ 


4> 


L 


CO 


Sur^cxce 

Figure 2. One-Sided Airfoil 
Equation (32) now becomes 


cW 




*X 


(1-33) 


The velocity potential perturbation can be related to the 
pressure potential perturbation by the equation 

j yp 

ijj [y ) eyp^cu>(A -^.)/\/^ cl A d-34) 

-oO 


Therefore, Eq.(33) can be replaced by an integral equation 
in the pressure potential perturbation l|j , which can be 
written as 





(1-35) 
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where both paths of integration are parallel to the ^ -axis, 

is the dummy variable for X , and is a point on this 
line whose X -coordinate is ■X . After rearranging the 
limits of integration, the area of integration is 

restricted to the lifting surface only, as shown in Figure 2, 
because a pressure differential cannot be sustained by the 
wake. The net pressure Ap equals — . 


The Green s Function and the Source Solution in the Uniform Case 


Under the conditions of uniform flow over a thin airfoil, 
the operators reduce to the form 


I- 

2 

V- 

C^- 

* vf 

C_ 

(1-36) 

X- 

2. 

V- 



(1-37) 

clearly 


X 


(1-38) 


Therefore JlL is Hermetian, also £ 
roles when the flow direction reverses. 

It can be shown that 




and 


£ 


exchange 


G r) - G*Tr,jp;> 


therefore , 

A ( jP> T)} = / r (G jO] ^(ir, jp3] 


(1-40) 
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and Eq. (22) can be reexpressed as 

<^V i c> ( 

giving, directly, the form of Green's function used in 
Eqs . (33) or (35) . 

However, itS(lf^jp) represents a pulsating source located 
at the point jp , it satisfies the equation 

=/ r Is (r,jp) & (r- jp^ d-42) 

thus ^ and Qr are identical. This fact has been assumed 
a-priori by most investigators, whereas, properly, it should 
have been derived from the Hermetian property of the operator 
£ in the uniform flow case. 


Conclusions 

Following the approach taken in the uniform flow case, it 
might appear correct to obtain, or at least to approximate, 
a source type solution, and to use it in place of the Green's 
function. However, since this would satisfy the equation 


^ XT i ^ “ S( fP-Jp) 


(1-43) 


it can be seen, by comparison with Eq. (30) that ^ and (j- 
different. Not only is the operator £ given in Eq. (20) 
different in form to the operator £ given in Eq.(31), but 
as described by Eq. (43) represents a source centered at 


are 
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the point Jp whereas as described by Eq. (30) represents 
a disturbance centered at the point IP . 

A physical explanation might be as follows. The cause of 
the potential perturbation is the normal flow or "downwash" 
perturbation on the surface. The result is the potential 
discontinuity. Therefore, it is incorrect to replace the 
potential discontinuity by a distribution of sources (in doublet 
form) because this reverses cause and effect. 

In conclusion, a ray tracing method based on using the wave 
equation, and on tracing from the "potential perturbation" 
point towards the "downwash" point, is bound to lead to 
incorrect results in the non-uniform flow case. However, 
if the adjoint equation is used, with the ray originating at 
the "downwash" point, and terminating at the "potential 
perturbation" point, realistic results might be obtained. This 
reverses the normally accepted roles of "disturbing" and 
"receiving" points. 



PART II 


APPLICATION OF MONTE CARLO METHOD TO 
ACOUSTICAL ANALYSIS 


II-l DISCUSSION 

Although analytical solutions to many acoustical 
problems have been known for many years, numerical tech- 
niques for the solution of acoustical problems involving 
arbitrary boundary conditions have lagged numerical develop- 
ments in other fields, for example, finite element methods 
in structural analysis. In fact one method of acoustical 
analysis which has shown some promise has been the develop- 
ment of finite fluid element methods, such as those by 

18 19 

Ziekiewicz and Chang and Oden 

An attempt was made, under this grant, to develop a 
Monte Carlo method for the solution of acoustical field 
problems involving arbitrary distributions of sources and 
boundary conditions. 

As a preliminary step, the proposed method was applied 3 
to a relatively simple problem involving a single sound 
source in a rectangular room. A reasonable degree of 
success was achieved, but it was found that many computer 
hours would be required to obtain accurate results. 

4 

The method was then applied to a problem involving a 
circular duct, for which an analytical solution was also 
available. Again, a reasonable degree of success was 
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achieved, but indications were that much more computer time 
would be needed for accurate results. 

In evaluating such results, it must be borne in mind 
that the Monte Carlo method is capable of solving problems 
involving arbitrary boundaries, for which no other solutions 
are available. Therefore, the need for further development 
depends on the need for such solutions. 

Two further lines of development remained unexplored 
at the termination of this work. These were 

(a) Adaptation of the method to moving media, such 
as within the duct or nozzle of a jet engine. 

(b) Adaptation of the method to problems involving 
diffractive scattering. This is a severe limitation, for 
example, although a rectangular room and a circular duct 
have been analyzed, a room with an open window, or an open 
ended duct, cannot be handled. 

The two applications mentioned above are discussed 
in more detail in the following sections. Listings of 
computer programs are to be found in Reference 4. 
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II-2 CONTEMPORARY WORK IN THE FIELD 

The first application of the Monte Carlo method to 
acoustics appears to have been by Allred and Newhouse 2 ^* 
who applied it to the determination of mean free paths of 
acoustical rays in rooms in order to estimate reverberation 
times. Their work contained an apparent procedural error 

22 9 ^ 

pointed out by Hunt , and by Schroeder , Schroeder and 
Kuttruff 2 ^J 2 ~*, applied the Monte Carlo method to the 
determination of the rates of occurrence of maxima in the 
response of a room. In a later work, Schroeder 26 showed 
that the standard deviation of the fluctuation obtained by 
this method agrees closely with analytical predictions. 

In a paper on the digital simulation of reverberations in 

27 

rooms, Schroeder proposes the use of the Monte Carlo 
method for tracing rays, with the suggestion that the 
amplitude of these rays be attenuated after successive 
collisions with the walls. Whereas, in the above references, 
the Monte Carlo method was employed to obtain statistical 
information on mean free paths of rays, and to estimate the 
occurrence of maxima, the present Monte Carlo application 28 
is directly aimed at obtaining solutions to acoustical field 
problems. Thus the end result is a time history of the 
pressure at a given point in the field. 
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II-3 APPLICATION TO RECTANGULAR ROOM PROBLEMS 
Introduction 

In attempting to set up a Monte Carlo method for 
acoustics, the author, with Thanedar, selected the direct 
representation of the acoustical pressure as the quantity 
to be determined. This does satisfy the wave equation, but 
is not a quantity which satisfies conservation laws. 

Therefore, for example, the pressure in a region is not 
directly proportional to the number of rays passing through 
it. However, an algorithm has been developed for local 
pressure which does depend on the characteristics of the 
rays passing through a small region, which is referred to as 
a "test cell". 

The process of ray tracing used in the Monte Carlo 
method starts with selection of a ray tube pulse of appropriate 
weight, originating at a source. Each ray, in turn, is 
traced through all of its interactions with the boundaries, 
which may result in reflection, or, if the boundaries are 
assumed to be absorbent, in the chance of a sudden death. 
Tracing of a particular ray is terminated when the required 
solution time has elapsed if it has not already suffered 
sudden death. During all of this time, the pressure-time 
history is accumulated at 'receiving' points by using the 
appropriate algorithm. 



Definitions 


C = Velocity of sound 


D{ 3- 
E{ 1- 

L 

I = 

* 

If = 

T - 

k = 

K = 
M = 

W = 

N - 
= 

P " 
Po = 

Q - 

R. = 
S - 

t - 
T P - 
V = 
V/ = 
A* = 

9 = 

r- = 
c r = 

-a = 


Variance of a quantity 

Expected value of a quantity 

Index on source pulse 
% 

Upper limit on U 
Index on time step 
Upper limit on ^ 

Index on sample 
Upper limit on k. 

Mass of air 

Index for ray count 

Upper limit on ft. 

Total number of pulses 

Total atmospheric pressure 

Static pressure 

Source strength 

Distance travelled by ray 

Index 

Time 

Time duration of a pulse 
Volume 

Statistical weight 
Increment of a quantity 
Density 

Time related to a pulse 
Standard deviation 
Solid angle 
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Calculation of Acoustical Pressures 

The rectangular room enclosure, with a single monopole 
source, is shown in Figure 3. A typical source time history 
of duration ) p is shown in Figure 4, this is approx- 
imated by square pulses of width A t , and of strength 
Q. L ' ^ = 1 to I , originating at time ( — C i— *fz~\ AT) 
It may be noted that Q has the dimensions of rate of 
change of volume flow rate, so that the volume flow rate 
will only return to zero if the integral of Q returns to 
zero, i. e. , if 

L f ( Qi- =f a Q(tUt = o (IlTl) 


where X - . Otherwise the mean pressure in an 

enclosure will steadily rise as air is injected into it* 

In the absence of boundaries, the acoustical pressure 
increment at a receiving point is given by 

Ap(t)- P (b) - Po = Q( t- R/c )/knR_ ,11-2) 

where Fv is the distance from the source to the point, 

POO is the total pressure, the ambient pressure, 
the velocity of sound and p ambient air density. When the 
boundaries form a rectangle, the effects of multiple 
reflections can be handled with relative simplicity, if 
the infiinte set r S — I to oO , represents all 
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of the possible distances travelled by an acoustical ray 
from the source to the point in question/ then the pressure 
at that point is 


co 

5 ~ I 


This result can be visualized more clearly by the 
method of images. If each wall, and each image of a wall, 
is replaced by its image, an infinite array of image rooms 
can be constructed, each containing a source, as shown in 
Figure 5. Then the pressure at the receiving point is the 
sum of the effects of all of these image sources, as well 
as the true source. 

Consider, now, the influence of some randomly selected 

ray pulse, let this be the pulse selected, and let the 

pulse have a strength and initial time which has 

been selected at random out of the values 'T ■ , l =• 1 to T. 

Further, let its influence be confined to a ray tube of solid 

angle , with randomly selected direction cosines, using, 

29 

for example, the method of selection given by the author 

(This method selects sets of direction cosines at random 
out of a population having a uniform probability density 
distribution with respect to solid angle. Thus spherical 
symmetry is preserved). 
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Fig. 5 Array of Image Rooms 
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Then, if the receiving point is contained within the 
ray tube for the time at a distance from the 

source, measured through all of the reflections, the effect 
on the pressure increment at the receiving point is represented 
by a S ” function, resulting in 

^phs(t) = AT 8(t7-%-^hs/cyifn-lT , hs (n-4) 


If time t is divided up into U” time cells of duration 

At • a nd if Ap >* is the mean value of APn.s(.b) taken over 

• ft w 1 

the (, ™ time interval, then, from Eq. 4 


A 




A Pnsft)it 

Ci-OA t 


- yArQn/^TrAt T? h5 

= O otherwise 


if the ray passes for the 
time during the ^ th interval 

(II-5) 


Although the method just described is under consideration 
for the case of curved boundaries, a simpler alternate method 
is possible for rectangular boundaries, in which a small 
fixed test cell of volume AV is set up, and is 

determined, the distance by which the ray penetrates 

the volume AV in the time interval. At the same time, 

the value of ' the mean distance from the source 

during penetration, is obtained. Then "R^ s ^ 

The two methods are shown in Figure 6. Clearly, AR, must 
be normalized by dividing by a characteristic length which 



RAY TUBE LOCATIONS 



PENETRATION 



TIME VOLUME AVAt 


Fig. 6 Alternate Methods of Accumulating 
Acoustical Pressure Contributions 




includes the ratio of relative cross section areas of the 
volume AV and of the ray tube, i.e., by 

= A.V / A-fL (ii-6) 

so that the pressure increment given by Eq. 5 is now 
written as 

A P ^ ^ p AT fan At AV <n-7) 

Finally, the expected pressure is obtained by averaging 
over all of the h4 pulses followed, where H = l to |\| , and 
by normalizing by multiplying by the total number of pulses 
H-p which could have been selected, i.e., by 

lx* l/llSl = hT (I1 _ 8) 

Thus the algorithm for the pressure increment in the 
^ time interval is obtained by averaging the pressure 
disturbances and then substituting from Eqs. 7 and 8, which 
leads to 

A P i = 2 ^p n \ 

u t\=\\ Q 

N 
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It may sometimes be desirable to present this algorithm 
in the dimensionless form 

M 

Ap./^c - - (Tp/Nc ZWM~)^Z.Q h &R h( j I 11 ' 10 * 

Calculation of Standard Deviations 

Since the results of a Monte Carlo calculation are 
statistical samples, it is important that we know their 
standard deviation. In order to do this, the calculation 
implied in Eq. 9 is made in K blocks, each consisting 
of Kj^ selected ray pulses. Each selection is given equal 
statistical weight, so that the statistical weight of the 
fc* batch is: 

Wfc = Nl £ /iSj (11-11) 

Now we modify Eq.9 slightly, and write 

% 

V A P (5>’ T p/^^V'A.r)^ZL(5nRv, ( j m- 12 ) 

It is clear that 

K 

= X (11-13) 

II - < 
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where 0^ is the standard deviation of the predicted 
acoustical pressure. 

A flow diagram showing the computer logic necessary 
to perform these calculations is given in Figure 7. 

Verification of Method Applied to Rectangular Room 

It is prudent to check the consistency of a method 
such as this by considering the results which would be 
obtained by solving very simple problems. The first solution 
considered was that of a source in an unbounded medium. 

The second was that of the ultimate average pressure to be 
expected in a rectangular room. Both are treated in the 
following sectiohs by evaluating the algorithm given in 
Eq.9 and then comparing with the known solutions. 


Source in Unbounded Medium 

The problem is illustrated in Figure 8. A source which 
has a constant strength Q over a total period Tp 
is surrounded by an annular test volume of mean radius and 

thickness Clearly, every selected ray penetrates 

this volume for a distance &^?y^=£^at a mean distance of 

= R. . The number of rays penetrating the control 
volume during the period is , and the 

test cell volume is 4 -tt R & R. . Substituting these values 
into Eq.9, the pressure increment becomes 



j> Q / 4^ 


(11-15) 




Fig. 7 Flow Diagram of Monte Carlo. Method 
Applied to Rectangular Room Problem 
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Fig. 8 Source in Unbounded Medium 
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which is the known analytical solution for the problem. 


Source in Bounded Medium 

The total air mass introduced into the room by the 


source is equal to 


t r T 


(11-16) 


**0r) = sJ o J 0 Q(^)^dr 

= j^rJo QtpO'VLif X r Qlr)dr 

If the source strength Q(t) meets the condition given 
by Eq.l that it should cut off at time *T p , i.e., 

r t 

J 0 = O * t > T p (11-17) 

then, for all times the final mass introduced into 

the room becomes 

r T P 

AM = - y J TQ(t)cL T" 


(11-18) 


and, provided that this is small compared to the mass O'Y 
already in the room, where \/ is its volume, the expected 
value of the adiabatic pressure rise in the room, E{Ap(t)} 
is equal to 


Cb) 1 = c :V m/V 


= - (j>e/v) J ‘T'Q(T)dT 3 t>Tp 


(11-19) 
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If is described in terms of equivalent pulses, 

as in Figure 4, then for t”*> If 

6 

E&p(tp] = -fj>cAT/\/).Z Ti.Qi. di-20 

L'— i 

Now consider the expected adiabatic rise predicted by 
the Monte Carlo method/ using Eg. 9. Because selection of 
any one of the £ pulses is equally possible, it is 
readily shown that 


E } = (fr, f N m Lb ) E [ Z ( Q n Af? h ^ ] 

I 

= Cj)c£^T/AVAb) E£ar} Ziti -rJQ t (11-21) 

The expected value of the penetration distance E^AR^is 
best obtained by using the method of images as shown in 
Figure 9, and by considering the effect of all of the image 
sources. Th’e total number of such images which can influence 
Ai is the number in an annular spherical volume of 

radius C Lt* ^ — ^nd thickness CAt • Since there is 
one source per room of volume V , the total number is 

3 ^ 

4rr c (tj -Tj At /V 


This number is actually equal to the number of 
reflections in the room during the corresponding period. 




RECEIVING POINT 


Fig. 9 Average Pressure by Image Method 
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The expected penetration into the test cell &V of any 
ray from one of these image sources 

(aV I ( tj-rS&t)* ( cAt) 

Therefore, the expected value of the penetration 
per pulse is the product of the above with the number of reflections 

= CAV &£: /y (11-22) 

Substituting into Eq.21, and assuming the cutoff 
condition given in Eq.l j- 

E - } - (ycWv3 

X 

= - (pc at/v) Z T- l Q I (II - 23) 

1=1 

which is in agreement with Eq.20 obtained by considering 
adiabatic compression. 


Calculations for Rectangular Rooms 
Double Rectangular Pulse 

Two sets of calculations are presented here. The 
first assumed a double rectangular pulse of 0.2 seconds 
duration as shown in Figure 10. A drawing of the rectangular 
room, specifying the essential dimensions is shown in 
Figure 11. The supporting medium was considered to be air 
with a mass density equal to 0.002378 slugs/ft and 

a velocity of sound C equal to 1100 ft/sec. Taking the 
space-time volume (i.e.&V&fcT) Q f the test cell to be 
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3 

0.4259 ft -sec, five blocks of data were obtained including 
a total of 35000 selected ray pulses. These calculations 
are presented for the duration of the pulse in Figure 12, 
where results from the five blocks are shown separately, 
and the overall mean values are shown plus and minus one 
standard deviation. One would expect the correct result 
to be within these limits about 70% of the time. For 
comparison, the results obtained by Mintzer 30 using the 
Laplace transform approach are given in Figure 13. 

One block of calculations for 5000 ray pulses was 
continued out to 0.4 seconds, twice the pulse duration. These 
results are shown in Figure 14, where they are compared with 
the expected adiabatic pressure rise as predicted by Eq.23. 

One Cycle Sinusoidal Pulse 

The second set of calculations was prepared for a 
comparison with calculations by the normal mode method. In 
order to make the results more convergent, a sinusoidal 
pulse was assumed, as shown in Figure 15. The Monte Carlo 
results are shown compared with the normal mode calculations 
in Figure 16 . These calculations consisted of ten blocks 
totalling 50000 selected ray pulses, and took 211 seconds of 
core time on the CDC 6400 computer at the University of 
Virginia’s Computation Center. The normal mode solution 
included all modes from (1,0,0) up to 0.5,15,15) and required 
118 seconds. Therefore, computation times were comparable. 
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Fig. 12 Pressure Time History for Square 
Pulse in Rectangular Room 
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Fig. 16 Pressure Time History for Single Sinusoidal 
Pulse in a Rectangular Room 
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It will be seen that the results of the two methods 
are in general agreement, when it is considered that the 
normal mode solution should only be expected to be within 
the limits of plus or minus one standard deviation about 
70% of the time. 



55 


II-4 APPLICATION TO CIRCULAR DUCT PROBLEMS 

4 

A solution to the problem of a source at the closed 
end of a circular duct was obtained by a method which might 
be described essentially as follows. 

A central ray was selected as in the rectangular 
room problem, and was traced through collisions with the 
cylindrical walls of the duct. In addition, two rays were 
selected making small angles with the central ray, so that 
the intersection of the three rays with a normal plane formed 
the vertices of a right angled triangle. These additional 
rays were traced. 

The test volumes were small spheres, to simplify the 
calculation of the penetration distance Af?. . Whenever a 
test volume was penetrated, the intersections of the rays 
with a normal plane were found, and the area of the triangle 
so formed was calculated. The distance which would have 
been travelled by the originating bundle of rays to form 


the same area was calculated, and taken as R . The value 

so calculated was substituted in place of in Eq. (II-9) 

o 

A typical reflection of the three rays from a curved 
surface is illustrated in Figure 17, while a comparison of 
a Monte Carlo solution with a normal mode solution is shown 

in Figure 18. In both cases?, a source on the center of the 


closed end of a 7.5" radius duct was assumed to emit for 


one cycle at a period of .002 seconds. The pressure was 
obtained at a point 10" from the end, and 3.16" from the 


axis . 
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Fig. 17 Incident and Reflected Rays 
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Fig. 18 Pressure vs Time for a Single Cycle Source at 
end of a Cylindrical Duct (Comparison of the 
Monte Carlo and the Normal Mode Solutions) 
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In these calculations/ a total of 2500 ray pulses 
were followed, divided into five equal blocks for the 
purpose of estimating deviations. These took 189 seconds 
of core time on the CDC 6400, or ,076 seconds per pulse. 

In contrast, the rectangular room solution required 211 
seconds for 50000 ray pulses, or .0042 seconds per pulse. 
The nearly 20x1 time increase was due to the much greater 
complexity of handling the wall reflections, as well as the 
to the fact that three rays had to be traced in place of 


one. 
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